struct distribution_ggx_parameters
{
  vec2 m_Alpha;
  float m_Anisotropy;
  float m_AlphaPrime;
};

void distribution_ggx_initialize(out distribution_ggx_parameters oParameters, in float iRoughness, in float iAnisotropy)
{
  oParameters.m_Alpha.x = vPow2(iRoughness);
  oParameters.m_Alpha.y = vPow2(iRoughness * (1.0 - iAnisotropy));
  oParameters.m_Anisotropy = iAnisotropy;
}

void distribution_ggx_initialize_area(out distribution_ggx_parameters oParameters, in float iRoughness, in float iAnisotropy, in float iAlphaPrime)
{
  oParameters.m_Alpha.x = vPow2(iRoughness);
  oParameters.m_Alpha.y = vPow2(iRoughness * (1.0 - iAnisotropy));
  oParameters.m_Anisotropy = iAnisotropy;
  oParameters.m_AlphaPrime = iAlphaPrime;
}

float distribution_ggx_g(in distribution_ggx_parameters iParameters, in vec3 iN, in vec3 iV, in vec3 iL, in vec3 iX, in vec3 iY, float iNoL, float iNoV)
{
  // Vector projection: a_proj = a - (a*n)/(n*n) * n
  //const float length_n_sqr_inv = 1.0 / vDot(iN, iN);
  //const vec3 l_proj = iL - iNoL * length_n_sqr_inv * iN;
  //const vec3 v_proj = iV - iNoV * length_n_sqr_inv * iN;
  // l||n, dot(iL,iN)=1 and v||n, dot(iV,iN)=1 are special cases
  //const float cosphi_l_sqr = (iNoL > 0.9999) ? 0.0 : vPow2(vDot(iX, vNormalize(l_proj)));
  //const float cosphi_v_sqr = (iNoV > 0.9999) ? 0.0 : vPow2(vDot(iX, vNormalize(v_proj)));
  //const float alpha_x = iParameters.m_Alpha.x;
  //const float alpha_y = iParameters.m_Alpha.y;
  //const float a_l = cosphi_l_sqr * alpha_x + (1.0 - cosphi_l_sqr) * alpha_y;
  //const float a_v = cosphi_v_sqr * alpha_x + (1.0 - cosphi_v_sqr) * alpha_y;
  //const float s1 = vSqrt(vMax(0.0, 1.0 - a_l + a_l/vPow2(iNoL)));
  //const float s2 = vSqrt(vMax(0.0, 1.0 - a_v + a_v/vPow2(iNoV)));
  //return 2.0 / (s1 + s2);

  const float alpha_x = iParameters.m_Alpha.x;
  const float alpha_y = iParameters.m_Alpha.y;

  //drop the anisotropic version for now... too many issues

  //const float length_n_sqr_inv = 1.0 / vDot(iN, iN);
  
  //vec3 v_proj = iV - iNoV * length_n_sqr_inv * iN;
  //vec3 l_proj = iL - iNoL * length_n_sqr_inv * iN;
  //
  //const float v_proj_len_sqr = vMax(1e-3, vDot(v_proj, v_proj));
  //const float l_proj_len_sqr = vMax(1e-3, vDot(l_proj, l_proj));
  //
  //l_proj = l_proj / sqrt(l_proj_len_sqr);
  //v_proj = v_proj / sqrt(v_proj_len_sqr);

  //const float cosphi_v_sqr = vPow2(vDot(iX, v_proj));
  //const float cosphi_l_sqr = vPow2(vDot(iX, l_proj));

  const float a_v = alpha_x;//cosphi_v_sqr * alpha_x + (1.0 - cosphi_v_sqr) * alpha_y;
  const float a_l = alpha_x;//cosphi_l_sqr * alpha_x + (1.0 - cosphi_l_sqr) * alpha_y;

  const float s1 = vSqrt(vMax(1e-6, 1.0 - a_l + a_l/vPow2(iNoL)));
  const float s2 = vSqrt(vMax(1e-6, 1.0 - a_v + a_v/vPow2(iNoV)));
  
  return 2.0 / (s1 + s2);
}

float distribution_ggx_d(in distribution_ggx_parameters iParameters, in float iNoH, in vec3 iH, in vec3 iX, in vec3 iY)
{
  const float alpha_x = iParameters.m_Alpha.x;
  const float alpha_y = iParameters.m_Alpha.y;
  const float alpha_x2 = alpha_x * alpha_x;
  const float alpha_y2 = alpha_y * alpha_y;
  const float NoH2 = iNoH * iNoH;
  const float XoH2 = vPow2(vDot(iX, iH));
  const float YoH2 = vPow2(vDot(iY, iH));

  if (iParameters.m_Anisotropy > 0.0)
    return (1.0 / (PI * alpha_x * alpha_y * vPow2(XoH2/alpha_x2 + YoH2/alpha_y2 + NoH2)));

  return (alpha_x2 / (PI * vPow2(NoH2 * (alpha_x2 - 1.0) + 1.0)));
}

vec3 distribution_ggx(in distribution_ggx_parameters iParameters, in fresnel_dielectric_parameters iFresnelParameters,
                      in vec3 iN, in vec3 iV, in vec3 iL, in vec3 iH, in vec3 iX, in vec3 iY,
                      in float iNoL, in float iNoV, in float iNoH, in float iLoH)
{
  const float d = distribution_ggx_d(iParameters, iNoH, iH, iX, iY);
  const float g = distribution_ggx_g(iParameters, iN, iV, iL, iX, iY, iNoL, iNoV);
  const float m = (d * g) / vMax(1e-6, 4.0 * iNoL * iNoV);
  const vec3 f = fresnel_dielectric(iFresnelParameters, iLoH);

  return max(vec3(0.0), m * f);
}

//No-fresnel version
float distribution_ggx(in distribution_ggx_parameters iParameters,
                      in vec3 iN, in vec3 iV, in vec3 iL, in vec3 iH, in vec3 iX, in vec3 iY,
                      in float iNoL, in float iNoV, in float iNoH, in float iLoH)
{
  const float d = distribution_ggx_d(iParameters, iNoH, iH, iX, iY);
  const float g = distribution_ggx_g(iParameters, iN, iV, iL, iX, iY, iNoL, iNoV);
  const float m = (d * g) / vMax(1e-6, 4.0 * iNoL * iNoV);

  return max(0.0, m);
}

//Area light version
vec3 distribution_ggx_area(in distribution_ggx_parameters iParameters, in fresnel_dielectric_parameters iFresnelParameters,
  in vec3 iN, in vec3 iV, in vec3 iL, in vec3 iH, in vec3 iX, in vec3 iY,
  in float iNoL, in float iNoV, in float iNoH, in float iLoH)
{
  float d = distribution_ggx_d(iParameters, iNoH, iH, iX, iY);
  d = max(d * iParameters.m_AlphaPrime, 0.0);
  const float g = distribution_ggx_g(iParameters, iN, iV, iL, iX, iY, iNoL, iNoV);
  const float m = (d * g) / vMax(1e-6, 4.0 * iNoL * iNoV);
  const vec3 f = fresnel_dielectric(iFresnelParameters, iLoH);
  return max(vec3(0.0), m * f);
}

float distribution_ggx_pdf(in distribution_ggx_parameters iParameters, in fresnel_dielectric_parameters iFresnelParameters,
                          in vec3 iN, in vec3 iV, in vec3 iL, in vec3 iH, in vec3 iX, in vec3 iY,
                          in float iNoL, in float iNoV, in float iNoH, in float iLoH)
{
  const float d = distribution_ggx_d(iParameters, iNoH, iH, iX, iY);
  return d * iNoH / (4.0 * iLoH);
}

vec3 distribution_ggx_sample(in distribution_ggx_parameters iParameters, vec2 iRand)
{
  const float phi = iRand.x * 2.0 * PI;
  const float a = iParameters.m_Alpha.x; //Isotropic only
  const float value = (1.0-iRand.y)/(1.0+(a*a-1.0)*iRand.y);
  const float cosTheta = vSqrt(value);
  const float sinTheta = vSqrt(1.0 - value);
  const vec3 H = vec3(sinTheta * vCos(phi), sinTheta * vSin(phi), cosTheta);
  return H;
}
